%%Getting declination data from US historic calculator
% for the whole united states


% Bounding coordinates
min_lat = 26;
max_lat = 50;
min_long = -126
max_long = -64;
step_size = 0.5;
start_year = 1750;
end_year = 1985;
step_year = 1;

baseurl = 'http://www.ngdc.noaa.gov/geomag-web/calculators/calculateUshistoric?';

data = zeros([(max_lat-min_lat)/step_size, (max_long-min_long)/step_size,(end_year-start_year)/step_year]);
lat_dim = 1;
for latitude = min_lat:step_size:max_lat,
    lon_dim = 1;
    for longitude = min_long:step_size:max_long,
        year_dim = 1;
        for year = start_year:step_year:end_year,
       
        queryString = [baseurl 'lat1=' sprintf('%5.2f', latitude) ...
            '&lon1=' sprintf('%5.2f', longitude) '&startYear=' sprintf('%4d', year)...
            '&endYear=' sprintf('%4d', year) '&resultFormat=xml'];
        
        block = urlread(queryString);
        
       % while  any(strfind(block, 'The calculator seems to be down at the moment'))
             block = urlread(queryString);
       % end;
        pause(0.01);
        
        A = strfind(block,'<declination>');
        B = strfind(block,'</declination>');

        if (~isempty(A) & ~isempty(B) & length(A) == length(B)),
            
               data(lat_dim,lon_dim,year_dim) = str2num(block(A(1)+14:B(1)-1));
               
        else,
           % display(block);
        end;
        
        year_dim = year_dim + 1;
        end;
        
        lon_dim = lon_dim + 1;
        end;
        lat_dim = lat_dim + 1;
        display(latitude);
end;

%%
    
R = georasterref('RasterSize', size(squeeze(data(:,:,1))),'LatLim',[min_lat max_lat],'LonLim',[min_long max_long]);       
load coast;
year = start_year;
load /home/mnair/projects/geomag/colortable.mat;
for i = 1:48,
    
    grid2image(squeeze(data(:,:,i)),R);
    set(gcf, 'papersize', [6 3]);
    set(gcf, 'PaperPosition', [0 0 6 3]);
    colormap(colortable);
    set(gca,'FontSize',10);
    h=get(gca,'XLabel');
    set(h,'FontSize',10);
    h=get(gca,'YLabel');
    set(h,'FontSize',10);
    caxis([-30,30]);
    h=colorbar;
    set(h,'FontSize',10);
    title(sprintf('Magnetic Declination for Year %4d',year));
    hold on;
    plot(long,lat);
    saveas(gcf,['/home/mnair/projects/geomag/ushistoric/' sprintf('%4d',year)],'pdf');
    year = year + step_year;
    cla;
end;
   
%% Testing the programmatic access

%queryString = 'http://www.ngdc.noaa.gov/geomag-web/calculators/calculateDeclination?lon1=105&lon1Hemisphere=W&lat1=40.051051&lat1Hemisphere=N&startYear=2011&startMonth=11&startDay=14&resultFormat=xml';

queryString = 'http://wwwtest.ngdc.noaa.gov/geomag-web/calculators/calculateUshistoric?lat1=45&lon1=-120&startYear=1985&endYear=1985&resultFormat=xml';

nBadQry = 0;

for i = 1:1000,
    
block = urlread(queryString);

if any(strfind(block, 'The calculator seems to be down at the moment'))
             nBadQry = nBadQry + 1;
        end;
end;